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ABSTRACT 



A computer simulation of a correlator receiver was 
developed and exercised to study the impact of a model-based 
signal processing algorithm on the detection of transmitted 
CW and LFM pulse acoustic signals incident on a planar array 
of electroacoustic transducers. The model of the ocean 
communication channel incorporates a space-variant sound 
speed profile. The transducer output electrical signals are 
cophased by an FFT beamformer via phase weighting, and 
summed to form a total array output signal. The total array 
output signal is correlated with a delayed replica of the 
transmit waveform and compared to a Neyman-Pearson thres- 
hold. Receiver performance is measured using a Monte Carlo 
technique to estimate the probability of detection for a 
fixed probability of false alarm versus the signal-to-noise 
ratio at the input of a single transducer. White, zero-mean, 
Gaussian transducer noise is assumed to facilitate compari- 
son between theoretical and simulated performance. Results 
indicate that model-based signal processing provides signi- 
ficant improvement of receiver performance. 
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INTRODUCTION 



I . 

Model-based signal processing is described by Mendel 
[Ref. 1] as an approach that exploits knowledge of the 
underlying physics of a problem to develop signal processing 
algorithms. Use of the approach implies that some a priori 
knowledge exists regarding the problem under consideration. 

In the case of an underwater acoustic communication problem, 
such a model has been developed by Ziomek [Ref. 2]. 

Ziomek derived a time- invariant , space-variant transfer 
function of the ocean volume based on the WKB approximation, 
which is an approximate solution of the linear, inhomogeneous, 
scalar wave equation describing the propagation of small- 
amplitude acoustic pressure waves when the speed of sound is 
a function of depth. Based on the transfer function of the 
ocean volume, Ziomek [Ref. 2:pp. 257-261] also derived an 
equation describing the output electrical signal at each 
element of a planar array of point sources. The output 
signal is described in terms of the frequency spectrum of 
the transmitted electrical signal, the transmit and receive 
planar arrays, and the random ocean medium transfer function. 
Vos [Ref. 3] used these results to develop a computer program 
that generates time-samples of a real baseband output 
electrical signal at each element in a receive planar array 
as a function of a variable ocean medium sound speed profile. 
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planar array size, array far-field beam patterns and func- 
tional form of the transmit signal. Ziomek has since modi- 
fied this program to generate time-samples of the complex 
envelope of real bandpass output electrical signals. 

The research documented in this thesis has the following 
objectives : 

develop a computer simulation of a correlator 
receiver which processes the output electrical signals 
generated by the computer program developed by Ziomek 
and Vos; 

apply the concept of model-based signal processing 
to the development of the signal processing algorithm 
used by the receiver; 

determine the effectiveness of the approach in the 
detection of signals from a planar array of point 
source elements. 

Since the effects of the ocean medium on the signals 
processed by the receiver are embodied in the random ocean 
medium transfer function, the basic question to be addressed 
may be stated as follows. Can a priori knowledge of the 
ocean medium, based on physical principles of acoustic wave 
propagation, be used to improve the detection of signals in 
a receiver, and to what degree? 

Ziomek' s use of linear systems theory to develop a 
transfer function model of the ocean medium immediately 
suggests the use of a compensating filter at the array 
output to remove the undesirable time delays due to system 
geometry and wave propagation effects. This filter would 
ideally cophase the signals at each element in the planar 
array, resulting in maximum signal output when the signals 
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are added together. This filter is implemented in the 
frequency domain through the use of Discrete Fourier 
Transforms (DFT) , and the approach is exactly analogous to 
the FFT beamforming procedure discussed by Ziomek [Ref. 2: 
pp. 153-176] . 

In the frequency domain, the time delays due to system 
geometry and wave propagation effects are represented as 
phase shifts which may be eliminated if known a priori. The 
concept of model-based signal processing is applied here to 
obtain the proper compensating phase shift for the known 
system geometry and wave propagation conditions. 

Section II describes the theory used to develop the 
receiver model. The system context within which the 
receiver operates is described and related to previous 
investigations. A functional description of the receiver 
is shown, and each of the major functional blocks is explained 
in some detail. Finally, a statistical description of the 
receiver's performance is developed. 

The computer implementation of the receiver structure is 
described in Section III. The logical flow of the computer 
program is discussed and related to the receiver descrip- 
tion. The use of multiple trials to estimate the probability 
of detection is explained. Each of the major subprograms 
is characterized in terms of function and implementation. 
Verification of the computer simulation is discussed last. 
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Section IV presents the data obtained from the simula- 
tion when a rectangular-envelope, continuous wave (CW) 
pulse or a rectangular-envelope, linear-frequency-modulated, 
(LFM) carrier is transmitted. Receiver performance is 
described by plotting the probability of detection (Pd) , 
for a given probability of false alarm (Pfa) , as a function 
of the input signal-to-noise ratio (SNR) at each element 
in the receive array. Plots are provided for different values 
of Pfa, and show the relative improvement in receiver per- 
formance as various medium and wave propagation effects 
are compensated for by the model-based, signal processing 
algorithm. In each plot, the receiver performance predicted 
by theory when all array element output signals are precisely 
cophased, and the array element .input noise is zero-mean, 
uncorrelated and Gaussian is shown as a dashed line. The 
dashed line is plotted from data obtained from a closed 
form expression relating Pd to Pfa as a function of array 
element input SNR, and is superimposed on the output plots 
to provide a baseline for judging the validity of the receiver 
simulation output data. 

Conclusions and recommendations are discussed in 
Section V. 
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II. THEORY OF THE RECEIVER MODEL 



A receiver operates in the context of a total communi- 
cation system consisting of a signal source (transmitter), 
a signal propagation medium (channel) and a signal sink 
(receiver) . It is the model of the communication channel 
that is of initial interest since the signal processing 
algorithm will depend in large part on the physics des- 
cribing the propagation of the signal through the channel. 

A. OVERVIEW OF THE COMMUNICATION SYSTEM 

Ziomek's model (Ref. 2] of the ocean medium is described 
in general as a time-variant, space-variant, random filter 
(transfer function) in which the index of refraction, or 
equivalently, the speed of sound, is a function of depth, 
and includes both a deterministic and a random component. 
However, in describing the electrical output signals from 
the receive aperture, the model becomes more restrictive in 
the sense that the channel is considered to be time-invariant, 
but still space-variant. Furthermore, the transmit and 
receive apertures are taken to be rectangular, planar arrays 
whose elements consist of complex weighted point sources. 
Complex weighting of the array elements provides the means 
for amplitude shading and beam steering both the transmit 
and receive array patterns . The complex weights are ideal 
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for removing the undesired effects of the channel on the 
output electrical signals from the receive array elements, 
and become the tool for applying the model-based, signal 
processing concept. Figure 1 depicts the geometry of the 
transmit and receive arrays in the transmission medium. 




Figure 1. System Geometry 



Ziomek goes further than simply considering system 
geometry in generating the output signals at the receive 
array. The spatial variance of the transfer function per- 
mits modeling of propagation effects due to the dependence 
of the index of refraction on depth. Figure 2 shows the 
ray-bending effects of propagation through such an 
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Figure 2. Ray Path Bending Due to Inhomogeneous Medium 



inhomogeneous medium. The phase shifts due to system geometry 
and propagation effects may be considered together or 
separately in the generation of the output electrical signal 
data in the computer simulation due to Vos [Ref. 3]. It 
is this output electrical signal data which is used as the 
input signal to the receiver. 

B. FUNCTIONAL DESCRIPTION OF THE RECEIVER 

The output electrical signal data processed by the 
receiver are time samples of the complex envelope of the 
bandpass electrical signal at the output of each element in 
the receive array. This implies that the bandpass acoustic 
signal incident on the receive array has been converted to 
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an output electrical signal by each transducer element in 
the array. Each real, bandpass electrical signal then 
passes through a quadrature demodulator to become an 
equivalent baseband, complex envelope signal, which is 
time-sampled and converted from analog to digital form. The 
complex envelope is represented by the I-channel (in-phase) 
and Q-channel (quadrature-phase) components generated by the 
quadrature demodulator. Time-sampling is done in a manner 
that satisfies the Nyquist criterion for the baseband infor- 
mation contained in the I and Q channels. Thus, many of 
the components associated with a receiving system are already 
contained within the simulation that generates the planar 
array output signal data. The receiver simulation assumes 
these components exist, and essentially provides signal 
processing of the complex envelope of the output bandpass 
electrical signals. The major functional blocks of the 
receiver include: 

- an array signal processor, 

- a correlator implementation of a matched filter 
receiver, 

- a magnitude-square operation, and 

- a threshold decision operation. 

Figure 3 shows the functional block diagram for the receiver 
model. The magnitude-square and the threshold decision 
operations are considered to be part of the correlator 
(matched filter) detector function block shown in Figure 3. 
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PLANAR 

ARRAY 




Figure 3. Receiver Block Diagram 

The receive array is assumed to be a rectangular planar 
array of M *N elements where both M and N are odd numbers. 
Each element is assumed to be an omnidirectional point 
source. The geometry of the planar array and associated 
mathematical notation is shown in Figure 4. 




Figure 4. Planar Array Geometry 
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The function of each quadrature demodulator is to 
convert the amplitude and angle-modulated , bandpass, 
electrical signal at each transducer output into its 
baseband, complex envelope. Thus, the low-pass complex 
envelope may be sampled at a much lower rate. Note that the 
use of a quadrature demodulator does imply that the carrier 
frequency of the transmitted signal is known. It should be 
emphasized that the output electrical signal at each element 
in the receive array is passed through its own quadrature 
demodulator before array processing begins. The quadrature 
demodulator is shown schematically in Figure 5, where the 
complex envelope of y(t), denoted y(t), is given by. 



y (t) = y c (t) + j y s (t) 



( 2 . 1 ) 




A RRAY 

TRANSDUCER 

ELEMENT 




Figure 5. Array Element Quadrature Demodulator 
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where : 




is the I-channel (in-phase) component 
of y ( t ) , and 




is the Q-channel (quadrature-phase) 
component of y(t). 



The array signal processor, shown schematically in 
Figure 6, is a FFT beamformer. The function of the array 
processor is to maximize the total output signal when the 
signals from each element in the array are summed. 



Figure 6. Array Processor Block Diagram 

The array processor is essentially a filter that compen- 
sates for system geometry and wave propagation effects on 
the signal transmitted through the channel. The array 
processor is implemented in the frequency domain where 




y'(lT s ,md x ,nd y ) 7(qf Q ,md x) nd y ) 



COMPLEX 

WEIGHT 

GENERATOR 
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filtering can be obtained by multiplication of each spectral 
component of the complex envelope by a complex weighting 
coefficient. The complex coefficients are computed from 
a knowledge of the channel model, thus the application of 
model-based signal processing. After applying the complex 
weights, the Inverse DFT is computed to recover cophased, 
time domain, complex envelope signals from each array ele- 
ment. By summing these cophased signals over all M x N array 
elements, a constructive interference effect is achieved, 
and the total output signal is maximized. 

The matched filter portion of the receiver is imple- 
mented by correlating the total time-sampled output signal 
and noise from the array, r(&T s ). Where, 

rUT ) = y T UT ) + n T UT s ) (2.2) 

with a time and frequency shifted replica of the complex 
envelope of the transmitted waveform, x(t). Since the 
phase of the received signal is, in general, unknown, the 
magnitude-square of the correlator output is taken as the 
input to the threshold detector. This input is compared to 
a preset threshold level y to determine the presence of a 
signal. The preset threshold y is computed from a Neyman- 
Pearson criterion. The schematic of the correlator/matched 
filter detector is shown in Figure 7, and has been shown by 
Van Trees [Ref. 4:pp. 244-247] to be the optimum receiver 
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7*it s ) 



Figure 7. Correlator/Matched Filter Detector 

for the detection of a bandpass signal, with random ampli- 
tude and phase, in the presence of white, Gaussian noise. 

C. STATISTICAL DESCRIPTION OF THE RECEIVER 

The output electrical signal data embodies several 
assumptions and restrictions that must be restated before 
further development of the receiver model. The signal data 
is assumed to be generated by a real, bandpass, acoustic 
field, y (t,r), incident upon an array of electroacoustic 

Pi — 

transducers. The acoustic wave, y M (t,r), is propagating 

/N 

in the in Q direction with velocity c, or 

y M (t,r) = x(t + [r-n o ]/c) (2.3) 



where ; 



and 



AAA 

r = xx + yy + zz 



( 2 . 4 ) 



✓\ AAA 

n = ux + vy+wz 
o o o o 



( 2 . 5 ) 



u = sin 0 cos Uj 
o o r o 



( 2 . 6 ) 



= sin 0^ sin ^ 
o o o 



( 2 . 7 ) 



w = cos 0 
o o 



( 2 . 8 ) 



Note that x(t) may be an arbitrary function of time. 

1 . Array Element Output Signal Description 

The acoustic field, incident upon an element in the 

array, is converted to an output electrical signal by the 

transducer, and is transformed into a baseband complex 

envelope at the output of the quadrature demodulator. After 

time sampling at a rate f = 1/T (in samples per second) , 

s s 

spatial sampling over the receive planar aperture in incre- 
ments of d and d , and assuming linearity in the transducer 
x y 

operation, the complex envelope of the output electrical 

th th 

signal at the m ,n element of the planar array lying in 
the XY-plane may be written as. 



y(2.T ,md ,nd ) 
J s x y 



= xUT 



s + 



[u o md x 



v 0 n< y /c] 



( 2 . 9 ! 
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For example, if the real, bandpass output electrical 
signal y(t) is an amplitude and angle modulated cosine wave 
of the form 



y(t) = a ( t) cos [2Tif c t + 0 (t) ] 



( 2 . 10 ) 



after quadrature demodulation, the complex envelope of y(t) 
may be represented as 



Y(t) = y c (t) + jy s (t) 



( 2 . 11 ) 



where : 



y ( t ) = a ( t) cos 0 ( t ) 

c 



( 2 . 12 ) 



and 



y s ( t) = a (t) sin 0 ( t) 



(2.13) 



or, in magnitude-phase form 



a ( t ) 






+ y s (t> 



(2..14) 



and 



0(t) = Tan 1 [y s (t) /y c (t) ] (2.15) 
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It is also possible to show that the envelope function, 

E(t), of the real, bandpass signal y(t) is 

E(t) = | a ( t) | (2.16) 

The main advantage of working with the complex 
envelope form of the real bandpass signal is that the 
baseband complex envelope may be sampled at a rate f (in 
samples per second) determined by the bandwidth of the 
baseband modulating waveform, independent of the carrier 
frequency . 

It is further assumed that the complex envelope of 
the transmitted waveform x(t) can be represented exactly, 
over an interval of T q seconds, by a finite, complex, 

Fourier series such as. 



K 




(2.17) 



where the complex Fourier coefficient c can be written as. 



c 



q 



a exp [+j 0 ] 
q F J q 



(2.18) 



where: 



f Q : is the fundamental frequency in Hertz of 

the signal x(t) with period T = 1/f 
I T o o 

seconds, and 
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K: is the maximum number^of harmonics used to 

represent the signal x(t) . 



The complex Fourier coefficients c can then be 

q 

determined directly from the DFT with respect to time of 
x ( £T ) , or 



L' 



- T. I 



£=-L ' 



x(£T g )W^ £ 



(2.19; 



where : 



( 2 . 20 ) 



( 2 . 21 ) 



and : 



L: is the total number of time samples taken 

during the time interval T Q seconds where 
L is a non-negative, odd integer, 

T Q : is the fundamental period or data record 

length in seconds, and 

T : is the sampling period in seconds (note 

f s - 1/T S ) • 



W L = exp [ + j 2 tt/L] 
L' = (L-l)/2 



To satisfy the Nyquist sampling theorem, it can be 
shown [Ref. 2:pp. 164-165] that. 



L > 2K + 1 



( 2 . 22 ) 
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and , 



T < T /L 
s — o 



(2.23) 



2 . Generation of Complex Weight Phase Factors 

With these assumptions it can be shown [Ref. 2:pp. 
165-166] that the normalized DFT with respect to time of the 
complex envelope of the output electrical signal y(£T ,md , 
ndy) can be expressed as. 



Y SN (qf o' md x' nc y = c q exp[ + je m ]exp[ ? j27Tfu o m dx /c] 



•x exp [ + j<|) n ]exp[:f j27Tfv o ndy/c] (2.24) 



where : 



f = f c + qf ; -K , . . . , q , . . . , K (2.25) 



and : 



f : is the carrier frequency, 

c : is the constant speed of sound in homogeneous 

ocean medium at the receive array, 

u q : is the direction cosine of the wave propa- 

gation vector along the x-axis , equation 

(2.6) , and 

v q : is the direction cosine of the wave propaga- 

tion vector along the y-axis . See equation 

(2.7) . 
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Phase factors 0^ and <J> are due to the separable 

complex weight, c . That is, 

mn 

c mn = c m d n = a m ex P M 6 J b n ex P (2 ' 26) 

Notice that by a proper choice of 0 and in equation 
(2.24) where , 

0 = ± j 2 7T f u md /c (2.27) 

m J o x 



and 



(t> = ±j2iTfv nd /c 

n J o y 



(2.28) 



the phase shifts due to system geometry may be completely 
cancelled leaving only the complex Fourier coefficients c^ 
at each element in the planar array. It should be noted 
that the phase correction factors are functions of both the 
array element position (mdx or ndy) , and the frequency f of 
each spectral component of the input bandpass signal spec- 
trum. That is, the frequency is given by equation (2.25). 

Once the phase shifts due to geometry are eliminated, 
taking the inverse DFT with respect to frequency of the 
spectrum of the output electrical signal Y (qf q , md^ , nd^) 
will yield cophased time signals at each element. Summing 
all the cophased signals will result in a maximum total 



29 



output signal y T (£T s ) from the array processor. That 
is. 






M' 

l 

m=-M ' 



N 1 

I 

n=-N ' 



y(£T ,md ,nd ) 
2 s x y 



(2.29) 



where : 



M' = (M-l)/2 



(2 . 30) 



and 



N' = (N-l)/2 (2.31) 

This approach is extended to compensate for the 
deterministic signal phase shifts caused by transmission 
through an inhomogeneous ocean medium. Since the variation 
in the speed of sound c(y) is assumed to be a function of y 
(depth) alone, the calculation involves only the <{> phase 
factor of the complex weight. Ziomek''' has shown that for a 
sound-speed profile c(y) with a constant gradient g, a 
closed form expression for the deterministic component of 
phase shift is, 

k c 

*MD ■ + 2Tr { T [n D (y) ‘ 1J + iy n J <2 - 32) 

o r 



1 

[Ref . 



Extension of work by Ziomek based on expression for 0 
2 : pp . 263-268] . MD 
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where : 



n D ( y) = c 0 /c D (y) 



(2.33) 



C D (y) = c o + gAy n (2.34) 

Ay n = (y r -y o ) + nd y (2.35) 



k = 2irf /c 
o o 



(2 . 36) 



n D (y): is the space-variant (with depth only) 

index of refraction, 

c D (y) : is the speed of sound at depth y, 

c Q : is the speed of sound at the transmit 
array , 

k : is the wave propagation constant at 

the transmit array, 

g : is the gradient (slope) of the sound 

speed profile, 

y Q : is the depth of the center element of 

the transmit array, and 

y r : is the depth of the center element of the 

receive array. 



The negative of the deterministic medium phase 
factor, equation (2.32), is simply added to the system 
geometry phase factor, equation (2.28), to obtain the total 
complex weight phase factor in the y direction. 
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3 . Array Processor Output Signal Statistics 



Development of a receiver model also requires a 
specification of the noise environment in which the receiver 
operates. Any realistic model for the noise environment is 
extremely complex. However in order to form tractable 
theoretical results that can be reasonably approximated in 
a computer simulation, zero-mean, additive, white, Gaussian 
noise (AWGN) is assumed at the output of each array element. 
Remember that the array element output is also the input to 
the quadrature demodulator. 

The AWGN model permits derivation of a closed form 
expression relating Pd, Pfa and array input Signal-to- 
Noise Ratio (SNR) . This type of noise process can also be 
reasonably approximated by using computer generated pseudo- 
random number sequences from a standard Gaussian random 
number generator. By comparing the theoretically predicted 
receiver performance with the results of the simulation, 
verification of the computer implementation of the receiver 
can be achieved. Once verified, the computer simulation 
can test more realistic noise models with some confidence 
in the resulting data. 

The input SNR at a single element in the array is 
defined as. 



a E{ |y (£T ,md ,nd ) | 2 } 

SNR. = - ^ (2.37) 

E { | n ( £T , md , nd ) | ^ } 

o Y 
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where E{ } denotes the expected value of the quantity within 
the curly braces. 

If the random input signal y ( £T ,md r ,nd ) is ergodic, 

s x y 

then the mean-square value of the signal can be found by 
computing the time-average instead of the ensemble average,, 
that is, 

E{ |y (£T ,md ,nd ) | 2 } = < |y UT ,md ,nd | 2 > (2.38) 

o Y o .X. V 

where < ( • ) > denotes the time-average of the quantity within 
the parenthesis, or 



< ( * ) > 



T 

o 



V 2 



-T 0 /2 



C) at 



(2.39) 



The integral in equation (2.39) can be approximated 
for computer simulation purposes by. 




where the infinitesimal dt is approximated by, 

dt = T s 



(2.40) 



(2.41) 
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Because of the assumption regarding the use of a 



finite Fourier series to represent the original transmitted 
waveform, the input signal power at a single element in the 
receive array can be computed from the sum of the magnitude- 
square of the complex Fourier coefficients [Ref. 5:pp. 44- 
45] , that is , 



These coefficients are easily obtained for each array 
element by computing the DFT with respect to time of the 
complex envelope of the output electrical signal data at 
each element, or 



~ 2 

< I y ( £T ,md ,nd ) I > = Y |c I 

1 2 s x y 1 ' __ cr 

2 q=-K ^ 



2 



(2.42) 



q= 



c 



1 

L 




L ' 

I ( y UT s ,md x ,nd y )w! 



,qi 



L 



(2.43) 





a. 



mn 



2 



(2.44) 



Substituting equations (2.42) and (2.44) into equation (2.37) 
the expression for input SNR at element (m,n) may be written 



as , 



SNR. 

1 



K 

l 



q=-R 





(2.45) 



or, by rearranging the variance, equation (2.45) becomes 



°mn - I Ic | 2 /SN Ri (2.46) 

q=-K ^ 

The measured input signal power at an array element 
and a desired input SNR parameter value SNR^ can be used to 
obtain the noise power (variance) required to scale the 
output from the random number generator. The ability to set 
a desired input SNR value is necessary in order to test 
receiver performance. 

Because the DFT and IDFT are linear operations, the 
noise statistics at the output of the array processor are 
still Gaussian, and uncorrelated in both spatial and 
temporal coordinates. The total noise signal n T (^ T s ) at 
the output of the array processor may be written as, 

M ' N ' 

n (£T ) = l l n(£T ,md x ,nd ) (2.47) 

m=-M' n=-N' y 

2 

The variance of the total noise a T is equal to the mean 
square value of the total noise signal from the array 
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processor, and is also equal to the sum of the variances 
of the noise signals at each element since the noise process 
at one element is by assumption independent of the noise 
process at any other element in the array. The variance 
of the total noise signal is then 



E{ |n(£T ) | } 



M' 

l 



N' 

l 



m=-M ' n=-N 1 



mn 



(2.48) 



If the variance a „ at each element is the same for all 

mn 

elements in the array, or 



mn 



= a 



(2.49) 



then, the variance of the total noise can be written in 
terms of the variance of the noise at each element in the 
array by substituting equation (2.49) into equation (2.48) 
to obtain, 



2 2 
a T = MNa 



(2.50) 



or, in terms of mean square values. 



E{ |n T UT s ) | 2 } = MNE{ |n(£T s ) | 2 } 



(2.51) 
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2 

where E{|n(£T s ) | } is the noise power at an array element 

when the noise power is considered to be the same for all 
array elements (m,n). 

Both the output signal and output noise components 
at all elements in the array are baseband complex envelope 
signals with I and Q channel components. The bandwidth of 
the noise signal is set by the bandwidth of the low-pass 
filter in the quadrature demodulator to W Hertz. For the 
purpose of the simulation, W is always adjusted to include 
the highest harmonic of the complex envelope signal, that 
is, 

W = Kf (2.52) 

Knowing the bandwidth W and using the AWGN assumptions, 
the noise power spectral density N q at each element can be 
related to the variance of the noise process at each element 
by, 



N = oL/(2Kf ) (2.53) 

o mn o 

The SNR at the output of the array processor SNR^ is 
defined as the ratio of the total signal power to the total 
noise power, or 



SNR a 



A 



E{ |y T UT s ) I 2 } 
E{ |n T UT s ) | 2 } 



(2.54) 
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where , 



2 M' N' 2 

E{ | y (£T ) [ ) = E { | l l y(£T ,md ,nd )| } (2.55) 

1 s m=-M ' n=-N ' s x y 

For the case of perfectly cophased signals with identical 
signal power at each element, the right-hand side of equa- 
tion (2.55) reduces to 



E{ |y T UT s ) 


| 2 } = (MN) 2 E{ |y UT s ) | 2 } (2.56) 


where E{ 1 y (£T ) 
s 


2 

| } is equivalent to the time-average power 


at each element 


in the array, and E{|y T (£T s ) | } is equivalent 



to the time-average power of the total output signal from 
the array processor. Using the uncorrelated and equal 
array element input power noise assumptions, the expression 



for SNR can be 


rewritten by substituting equations (2.51) 


and ( 2 . 56 ) into 


equation (2.54), or 


snr a 


E { |yUT ) | 2 } 

= MN r = MNSNR. (2.57) 

E { |nUT g ) | } 1 


Since the array 


gain (AG) is defined as. 


AG = 


10 log 1Q [SNR A /SNR i ] (dB) (2.58) 
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the array gain in this case becomes, 



AG = 10 log Q [MN] (dB) (2.59) 

Note that for a 5 x5 element planar array, M xN = 25, and 
the array gain is 13.979 dB . 

4 . Hypothesis Testing and the Neyman-Pearson Criterion 
The correlator/matched filter detector portion of 
the receiver, Figure 7, is modeled as a binary hypothesis 
testing problem using the Neyman-Pearson decision criterion. 
The two hypotheses, and , are defined as, 

/ y T U T s ) + n T UT s ) : H ] _ 

r(£T s ) = < (2.60) 

V n UT ) : H 

is 0 



where, in terms of the transmitted waveform x(t), 



y T (aT s ) 



M' N' 

c J y x [ (£T -t ) ,md ,nd ] 

h,. L .. , s mn' ;.x y 

m=-M n=-N 1 



x exp [+j 27 r<j) A 5 ,T J 



(2.61) 



and n T (£T s ) is given by equation (2.47). If one assumes 
that all array element output signals are cophased and 
identical, equation (2.61) reduces to 
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(2.62) 



y T (£T g ) = cMNx (£T g -T A ) exp [ + j27T(j) A £T s ] 



where : 



c = a exp [ + j9] (2.63) 



and. 



x mn : is the actual time delay in seconds at 

each element in the array due to range 
separation between individual elements 
of the transmit and receive arrays , 

t a : is the actual time delay in seconds due 

to range separation between transmit 
and receive arrays when all signals are 
cophased , 

<j) : is the actual' doppler shift in Hertz due 

to relative motion between transmit and 
receive arrays, 

a: is the amplitude attenuation factor that 

is, in general, a random variable, 

0: is the generalized phase shift of the 

received signal with respect to the trans- 
mitted signal. (In general, it is also a 
random variable dependent on both spatial 
and temporal coordinates.), and 

Hq : is the null or noise only (no signal) 

hypothesis . 



Matched filtering is obtained by correlating the 
complex envelope of the total array output r(£T g ) with the 
complex conjugate of the processing waveform g(£T g ). The 
functional form of g ( £T ) can be written as, 

b 



~ ^ A /N 

g(£T ) = x ( £T -x ) exp [ + j 2 tt4)£T J 

s s s 



(2.64) 
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which is a time and frequency shifted replica of the 
transmitted signal x(t) where 

A 

T : is the estimate of time delay at the receiver, 

and 

A 

4* : is the estimate of the doppler shift at the 

receiver . 

For the purposes of this study, the estimates of 
time delay and doppler shift are assumed to be precisely 
correct, or 



T = T - T. = 0.0 

A 



(2.65) 



and 



4> = <f>^ - 4> = 0.0 



( 2 . 66 ) 



The correlation, or inner product t between two 
functions r(t) and g (t) , is given by 



+ 00 



l = <r(t),g(t)> = / r(t)g*(t)dt 



(2.67) 



which is approximated in the simulation using the trapezoidal 
rule approximation to the integral, that is. 



L’-l r UTjg* UTJ + rlU+l)T ]g*[U+l)T ) 

l = l ^ — — T (2.68) 

£=-L’ S 
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The magnitude-square of the correlator output is 
taken for two reasons. First, the phase of the carrier 
frequency waveform is, in general, unknown. As the phase 
of the received carrier varies with respect to the phase of 
the quadrature demodulator local oscillator (LO) signal, the 
output of the quadrature demodulator would vary from a 
maximum negative to a maximum positive value depending on 
the phase difference between the carrier and the LO which 
usually is taken to be a uniformly distributed random varia- 
ble between 0 and 2 tt radians. The change in polarity of 
the quadrature demodulator output would propagate through 
to the output of the integrator in the correlator/matched 
filter detector. Figure 7. Taking the magnitude of the 
integrator output ensures that the input to the threshold 
comparator will always be non-negative regardless of the 
phase difference between the carrier and LO waveforms. 

Second, when the array element noise statistics are assumed 
to be Gaussian, the square of the magnitude of the integrator 
output yields an input to the threshold comparator which 
can be described statistically by exponential density 
functions for both and signal hypotheses. As will 
be shown in the following derivations, the exponential den- 
sity functions can be easily integrated to obtain a closed 
form expression for Pd and Pfa in terms of the SNR at the 
input to the threshold comparator. The output of the 
magnitude-square operation is the sufficient statistic 
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on which the binary decision is made. That is, choose 

if. 



l 



(2.69) 



or choose if. 



(2.70) 



Y 



Assuming that the total noise n^,(t) is a baseband 

Gaussian process of bandwidth W, the conditional probability 

density functions (pdf's) of the magnitude-square of the 

correlator output with and without a signal present can be 

2 

shown to be exponential. The conditional pdf's are given 
by 



p ( | £ | 2 I H q ) = l/(2c 2 )exp[-|£| 2 /(2a 2 ) ] 



(2.71) 



and , 



p ( 1 1 



H i ) 



l/(2cr 2 )exp[-|£| 2 /(2a 2 ) ] 



(2.72) 



For zero mean, AWGN, and cophased, equal energy (power) 
signals at each array element output, the variances of 



2 

Derivations for the pdf's, expressions for Pd, Pfa and 
the decision threshold Y were provided by Prof. L. J. Ziomek 
in private communication. 
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the magnitude-square of the correlator output can be shown 
to be 



°o = E{ I £ N | 2} / 2 = MNN o E^/2 



(2.73) 



and 



°1 = E {| i s | 2 + K n | 2 }/2 



= [ (MN) 2 E{a 2 } 



X ( t , 4>) I + MNN E~]/2 
1 o x 



(2.74) 



where : 



t . 1 : 

JM 1 



N 



E{a 2 }: 



X ( x , <J> ) | : 



is the magnitude-square of the correla- 
tor output when the input to the 
receive array is noise alone, 

is the magnitude square of the correla- 
tor output when the input to receive 
array consists of signal alone, 

is the power spectral density level of 
the noise signal at each array element 
output , 

is the energy in the transmitted signal 
which for simulation purposes is defined 
to be equal to the energy E~ in the local 
processing waveform g(t), y 

is the mean-square value of the amplitude 
attentuation factor (note: E{a 2 } = a^ 

only when deterministic effects are 
considered) , and 



is the magnitude-square of the auto- 
ambiguity function. Note that |x((x,4>) 
E? | X (t , (f>) | 2 , or in our case where x = 
and <J> = 0, | X ( 0 , 0 ) | ^ = e 2 since 

I X N (0,0) I 2 = 1. [Ref. 2?pp. 190-191] 



2 



0 
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To calculate the threhold y for a desired Pfa and 
input SNR, the integral of the pdf for the null, or Hq 
hypothesis, is set equal to the desired Pfa, and the result- 
ing equation is solved for the threshold. The magnitude- 
square operation results in particularly simple threshold 
equation when the input noise is assumed to be white, zero 
mean, and Gaussian. Computing the Pfa yields 



+oo 



Pfa = / P ( I £ | 2 H n ) d | ^ | 2 



Y 



= exp l-y/ (2 o q ) ] 



= exp[-y/(MNN E~) ] 



(2.75) 



By solving equation (2.75) for the threshold y, 



y = MNNE- £n [1/Pfa] (2.76) 

O X 

Once the threshold is obtained, the correlator 
output pdf for the H-^ hypothesis may be evaluated to obtain 
the Pd. 



Pd 



f p(KI 2 

Y 



H 1 ) d | £ | 



2 



exp [-y/2a 



2 

1 



] 



exp{-y/ [ (MN) ^E{a}^ 



X (t , <f>) 



'+MNN E~ } J 
ox 



(2.77) 
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Since in our problem x = 0 and <J) = 0 , 



| X ( x , <J>) 



|X(0,0) 



= E- 



x 



(2.78) 



and if we let 



E{a }E- 



SNR | ^ | 2 = MN ( — 



x. 



(2.79) 



then , 



MNE { a }E- 



Pd = exp{-y/[MNN E~ ( — 

OX L\ 



+ 1 ) ] } 



(2.80) 



and 



Pd = exp {-y/ [MNN E~ ( SNR . ~ , 2 +1)]} (2.81) 

° x 111 

Substitution of equation (2.75) into equation (2.81) 
gives the desired closed form expression relating Pd, Pfa 
and the SNR of the magnitude-square of the correlator output, 
that is, 



1/ [1+SNR | j | 2] 
Pd = Pfa 11 



(2.82) 



This result agrees in form with the result given in Van 
Trees [Ref. 4:pp. 246-247] for a similar single channel 
receiver model. Figure 8 graphically depicts the relationship 
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1 



7 



Figure 8. Density Functions of the Magnitude- 
Sqaure Correlator Output 

between the conditional pdf's, the decision threshold y, the 
Pfa and the Pd. 

The correlator output SNR, equation (2.79), can be 
related to the input SNR at a single element in the array 
through the array gain, and a factor resulting from the 
slightly different definitions of input SNR and correlator 
output SNR. In terms of array element input signal energy 
, the input SNR given in equation (2.45) may be rewritten 



as , 



SNR. 

l 



( 1/T )E~ 

° y 

2 

a „ 



E { a 2 }E~ 

X 



(2.83) 
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where , 



= E{a }E ' 



(2.84) 



is the average received energy at a single element in the 
array due to the transmitted signal x(t), and from equation 
(2.53) 



2kN 



mn 



2kf N 
o o 



(2.85) 



The magnitude-square correlator output SNR can be written 
as , 



SNR | ^ | 2 



A E{ l^s 






( 2 . 86 ) 



where , 



E(U S | 2 } 



(MN) 2 E{a 2 }E~ 



(2.87) 



and. 



E{|£ n | 2 } = MNN o E- 



( 2 . 88 ) 



Substitution of equations (2.87) and (2.88) into equation 
(2.86) gives 
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E { a 2 }E 



SNR | ~ j 2 = MN 



x 



N 



(2.89) 



By substituting equation (2.85) into equation (2.83) and 
rearranging terms, equation (2.83) becomes 



SNRj^ 



E{a 2 }E~ 
x 

2KN 

o 



(2.90) 



By rearranging terms in equation (2.90) and substituting 
the result into equation (2.89), the desired relationship 
expressing the magnitude-square correlator output SNR in 
terms of the array element input SNR at a single element 
in the array is obtained. That is, 



SNR|2|2 = MN2KSNR ± (2.91) 

Writing both sides of equation (2.91) in dB form yields, 



SNR | £ | 2 (dB) 



10 log 10 (MN) + 10 log ( 2K) + SNR^ 
AG 



(dB) 

(2.92) 



where AG is the array gain given in equation (2.59). 

To summarize. Section II. C provides the equations 
needed to implement a computer simulation of the receiver 
structure described in Section II. B. The implementation of 
the computer simulation is the subject of Section III. 
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Ill . 



COMPUTER SIMULATION OF THE RECEIVER 



The computer program RCVR simulates receiver operation 
through straightforward application of the equations and the 
concepts developed in Section II. Written in FORTRAN, the 
computer program RCVR consists of a top level, or main 
program, and nine subprograms. The computer program will 
be explained from a functional viewpoint. That is, the 
algorithms used to implement the receiver simulation will be 
related to the theoretical development outlined in Section 
II, but translation of these algorithms into FORTRAN state- 
ments will not be discussed. The main program will be 
described first. The description of the main program will 
be followed by a detailed discussion of each subprogram. In 
addition to explaining the computer simulation, the methods 
used to validate the receiver simulation output data will 
be presented. 

A. TOP LEVEL PROGRAM DESCRIPTION 

The organization and logic flow of the top level, or 
main program, is shown in Figures 9a through 9e. The 
functions of the main program include: 

- initializing the simulation run-time environment, 

- invoking subprograms in the proper sequence to 
process the input signal data. 

- providing control logic and noise generation algorithms 
needed to measure receiver performance, and 
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READY 



SGNLGN 



AMPWGT 



Figure 9a. Program RCVR Flowchart 
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V 




PHSWGT 



ARYPRO 



Figure 9b. 



Program RCVR Flowchart 





Figure 9c. Program RCVR Flowchart 
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AWGN 



ARYPRO 



Figure 9d. Program RCVR Flowchart 
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Figure 9e. Program RCVR Flowchart 
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- selecting the form of the simulation output data. 
Initialization of the run-time environment involves reading 
input data from file to internal storage, generating the 
local processing waveform, measuring the time-average input 
signal power at each array element, computing the baseband 
signal bandwidth, and generating the complex weights used 
by the array processor. The input data is read by a call 
to the subprogram READY which returns a set of simulation 
parameters in COMMON storage and the complex envelope 
electrical signal data for each array element in matrix 
form. The local processing waveform is obtained by a call 
to the subprogram SGNLGN. The time-average signal power 
at each element (m,n) in the planar array is found by using . 
the fact that the original transmit signal was synthesized 
from a finite Fourier series, and applying equation (2.43) 
to obtain the complex Fourier coefficients, c . Once the 
complex Fourier coefficients at each element in the planar 
array are obtained, the time-average power is computed using 
equation (2.42). The baseband signal bandwidth W is com- 
puted using equation (2.52). The separable complex weights, 
c m and d n inequation (2.26) are generated in two stages. 
First, the real amplitude factors, a m and b n in equation 
(2.26), are computed by a call to subprogram AMPWGT. Second, 
the real phase factors, 0 m and $ n in equation (2.26) , are 
obtained by a call to subprogram PHSWGT. The complex weights 
are then computed by combining the amplitude and phase 
factors as shown in equation (2.26). 
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The first signal processing step is to generate the 
total array output signal in the absence of additive noise 
by a call to subprogram ARYPRO. Subprogram ARYPRO which 
uses the complex weights c m and d^ to cophase the planar 
array output signal returns a total array output signals. 

The total array output signal energy is computed using 
equation (2.40) for later use in calculating a SNR at the 
output of the array processor. The ratio of array output 
SNR to the input SNR defines the array gain, and provides a 
check on the validity of the data generated by the array 
processor algorithm. 

Receiver performance is measured by computing a relative 
frequency estimate of the Pd over the specified range of 
array element input SNR values when the Pfa is a known, 
time- invariant parameter. Using the Pfa and a range of SNR 
values specified by the programmer in a receiver control 
data block, the main program establishes nine input SNR 
values for which the estimate of the Pd will be computed, 
and determines the number of trials, or runs, to be used in 
computing the estimate. 

At point B in the flowchart of Figure 9c, the main 
program enters a loop that initializes the noise source, 
performs an array gain calculation, computes the detection 
threshold and then enters a second, inner loop where the 
correlation detection is done. After exiting the inner 
loop, the relative frequency estimate of the Pd is computed. 
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The outer loop increments through the input SNR values 
determined earlier, and exits the loop when the specified 
maximum input SNR value is reached. 

The first step in initializing the noise generator is 
to determine the noise variance (power) at each array element. 
The power level, or variance, of the noise samples for a 
particular element in the array is a function of the time- 
average signal power at each element and the array element 
input SNR. The variance can be found using equation (2.46). 

By controlling the power level of the noise process at 
each array element, the receiver performance can be measured 
over a range of array element input SNR values regardless 
of the time-average signal power level at a particular 
array element. 

The variance computed in this manner is the variance of 
the complex envelope baseband noise signal. The variance of 
the real baseband noise signal produced by the noise generator 
subprogram AWGN in the I or Q channel is one-half the vari- 
ance of the complex noise signal since the complex envelope 
is the sum of the two independent I and Q channel noise sig- 
nals. Therefore, the variance of the complex envelope of 
the noise is divided in half prior to being passed to the 
noise generating subprogram as a scale parameter. 

The variance of the total complex envelope noise signal 

2 

o T is obtained by summing all the noise variances of the 

2 

array element complex envelope signals, a , as indicated 
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by equation (2.48). A power spectral density of the total 

2 

noise signal is computed using equation (2.53) where o T is 

2 

substituted for a . 

mn 

The SNR at the array processor output is computed using 
equation (2.54) where the time-average power of the total 
array output signal is substituted for the mean square 
ensemble average, and the mean square ensemble average of 
the total noise is taken to be the same as the total noise 
signal variance. The array gain is computed using equation 
(2.58), and is held in internal storage for later output in 
tabular format. 

The detection threshold y is computed using equation 
(2.76) where the total noise power spectral density is 
substituted for the MNN q factor and the energy of the local 
processing waveform E~ is used instead of the transmit 
signal energy E~. For the purpose of this simulation however, 
the energy in the local processing waveform is the same as 
the energy in the transmit waveform. That is, E~ and E~ 
are equal. 

The second, or inner loop, of the main program begins 
at point D in Figure 9c. Within the inner loop, the complex 
envelope noise signals are generated for all array elements, 
the total noise output signal is generated by the array 
processor subprogram ARYPRO, the total signal and total 
noise are summed and correlated with the local processing 
waveform, and the signal detection decision is made. The 
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inner loop terminates when the number of iterations through 
the loop equals the number of trials allowed to form the 
estimate of Pd. 

The time-sampled, complex envelope array element noise 
data is obtained by repeated calls to the noise generator 
subprogram AWGN. Each call to AWGN returns a properly 
scaled pseudorandom number representing one sample of the 
noise process in the I or Q channel at a particular array 
element with L time samples taken for each of the I and Q 
channels. The noise signal data is stored in the same 
matrix form as the array element output electrical signal 
data. 

The complex envelope noise data is submitted to the 
array processor subprogram for processing in exactly the 
same manner as the input signal data. The result is an 
array total noise output signal. Processing the noise alone 
in this manner provides some gain in execution speed and 
provides the flexibility to estimate the probability of 
false alarm directly, if desired. The DFT and IDFT are 
linear operations, and the principle of superposition holds. 
Therefore, the addition of the total noise and total signal 
at the output of the array processor is equivalent to adding 
noise to the signal at each element prior to passing the 
data to the array processor subprogram. 

Correlation of the sum of the total signal and total 
noise with the local processing waveform is accomplished by 
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using the trapezoidal rule to approximate the integral of 
equation (2.67). The magnitude-square of the correlator 
output is obtained by taking the complex product 

\l\ 2 =1-1* (3.1) 

The threshold detection portion of the receiver is 
implemented by comparing the output of the magnitude-square 
operation to the decision threshold, using a simple IF-THEN- 
ELSE binary branch. The number of hits, or times the output 
exceeds the threshold, are counted and stored. The process 
of generating noise samples, and making a hit or miss deci- 
sion continues through a large number of trials. Since the 
correlation is done using the (signal plus noise) hypothe- 
sis, Pd can be directly estimated using the ratio 

A 

Pd = HITS/TRIALS (3.2) 

The first approach taken to determine the minimum number 
of trials required to estimate Pd was based on equation 
(2.82). After rearranging terms, equation (2.82) can be 
written as 



Pfa = Pd I1+SNR] (3 

The central idea was to compute an estimate of Pfa, using 
equation (3.3), from the relative frequency estimate of Pd 
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in equation (3.2). The algorithm would then terminate when 
the computed estimate of Pfa differed from the Pfa input 
parameter by some arbitrary small amount. However, the use 
of equation (3.3) was found to be a poor test for establish- 
ing when the algorithm should terminate, and would not 
terminate the algorithm for Pd values greater than about 
0 . 6 . 

A perturbation sensitivity analysis of equation (3.3) 
can be performed by computing the total differential of 
equation (3.3). That is, 



dPf a 



3Pf a 
3Pd 



dPd + 



9Pf a 
3SNR 



dSNR 



dPf a = [1+SNR] Pd SNR dPd + Pd 1 1+SNRj iln (Pd) dSNR (3.4) 



Assuming that the SNR is a constant, or equivalently, the 
differential dSNR is zero, dPfa becomes 



dPf a = [1+SNRJ Pd SNR dPd 



(3.5) 



or 



APf a = [1+SNR] Pd SNR APd 



(3.6) 



where 



dPf a 



APf a 



Pfa - Pfa 



(3.7) 



and 



dPd ^ APd = I Pd -Pd 



(3.8) 



Dividing both sides of equation (3.6) by equation (3.3) 
yields 



APf a 
Pfa 



[ 1+SNR] 



APd 

Pd 



(3.9) 



where APfa/Pfa and APd/Pd represent the fractional error 
between the actual and estimated values of Pfa and Pd, 
respectively. 

Such an analysis indicates that the percent error in the 
computed estimate of Pfa is linearly related to the percent 
error in the estimate of Pd. The constant of proportionality 
relating the error in the computed estimate of Pfa to the 
error in Pd estimate is equal to (1 + SNR) where the SNR is 
taken at the output of the magnitude-square operation. Thus 
for SNR values greater than approximately 1, or 0 dB, a small 
error in the relative frequency estimate of Pd is scaled to 
a larger error in the estimate of Pfa found by using equation 
(3.3). Note that for a 5 x5 element planar array and a 
bandwidth of 5 times the fundamental frequency as in the CW 
pulse case, the array element input SNR can be found by 
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equation (2.92), and is roughly 24 dB less than the SNR at 
the output of the magnitude-square operation. That is, a 
0 dB SNR at the output of the magnitude-square operation 
corresponds to an array element input SNR of -24 dB . The 
end result is an estimate of Pfa that diverges wildly from 
the specified Pfa parameter at SNR values over which the 
receiver operates. In addition, equation (3.3) was derived 
using certain assumptions regarding the statistics of the 
input noise that may not be precisely duplicated by the 
pseudorandom noise data generated by the computer program. 

If the assumptions regarding the use of zero mean, uncorre- 
lated, Gaussian noise are not satisfied by the pseudorandom 
noise source, equation (3.3) does not hold, and the algorithm 
would not be suitable for other input noise models. For 
these reasons, the use of equation (3.3) was abandoned in 
favor of an empirically determined fixed number of trials to 
estimate Pd. 

The number of trials needed to estimate Pd was found by 
assuming that the absolute minimum number of trials should 
be the reciprocal of the Pfa parameter. That is, if a Pfa 
of 0.01 is specified, at least 100 trials must be taken to 
allow at least one chance in one hundred of a false alarm 
occurring. Even though a relative frequency estimate of Pd 
is being computed, Pd is related implicitly to Pfa through 
the decision threshold, and the value of the Pfa parameter 
should be taken into account when attempting to fix the 
number of trials needed to estimate Pd. 
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The number of trials was empirically determined by 
running the receiver simulation program repeatedly with 
increasing multiples of the minimum number of trials, and 
observing the effect on the estimate of Pd. Depending upon 
the value of the Pfa parameter, it was found that four to 
eight times the minimum number of trials would produce 
curves that did not change appreciably as the number of 
trials was increased further. Therefore, the fixed number 
of trials used to estimate Pd was arbitrarily set at 10 x 1/Pfa 
for Pfa of 0.1, and to 5 x 1/Pfa for Pfa of 0.01. The smaller 
multiplier for the Pfa of 0.01 became necessary due to limits 
on computer resources. 

1 . Subprogram READY 

The function of subprogram READY is to obtain the 
simulation parameters and the complex envelope output 
electrical signal data from a data file. The data file is 
generated by the ocean communications channel simulation 
computer program developed by Vos and Ziomek [Ref. 3]. That 
is, READY forms the interface between the RCVR simulation 
and the ocean communications channel program. A flowchart 
of READY is shown in Figure 10. 

The simulation parameters are read in first, and 
are stored in the named COMMON blocks: HEADER, SIGNAL, 

ARRAY, and MEDIUM. The named COMMON blocks provide the 
mechanism for communicating simulation parameters into the 
various subprograms . The HEADER data documents the type of 
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Figure 10 . Subprogram READY Flowchart 

communication channel, and the date the data was generated 
by the ocean communication channel simulation. The SIGNAL 
data includes: the fundamental frequency f Q and the period 

T 0 ; the number of harmonics K; the same rate f and the 
sample period T s ; the number of time samples L; the carrier 
frequency f ; and the number of zeroes padded to the input 
signal data, if any. The ARRAY data includes: the number 

array elements M and element spacing d^ in the x-direction; 
the number of elements N and element spacing d in the y- 
direction; the direction cosines u r and v , representing the 
direction of the direct path from the transmit array to the 
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receive array; the direction cosines and v^, representing 

the direction the transmit array beam pattern was steered in 

the communication channel simulation; the depth y Q of the 

center element of the transmit array; and the depth y^ of 

the center element of the receive array. The MEDIUM data 

includes the speed of sound c Q at the center element of the 

transmit array; the sound-speed-profile gradient g; and 

the speed of sound c „ at the center element of the receive 

yr 

array . 

The time-sampled, complex envelope, output electri- 
cal signal data is read in next. The data is stored in the 
complex matrix variable YCE with dimensions L, M and 
N. The maximum values of L, M, and N are limited to 33, 11 
and 11, respectively. 

2 . Subprogram SGNLGN 

The function of subprogram SGNLGN is to generate the 
time samples of the complex envelope and to compute the 
energy of the local processing waveform given by equation 
(2.64). A flowchart of Subprogram SGNLGN is shown in Figure 
11 . 

The local processing waveform g(£T g ) is synthesized 
from a finite, complex Fourier series with provisions for 

A 

incorporating the estimates of time delay t and doppler 

A ^ 

shift <j> in the total received signal y(£x ). Accurate esti- 

b 

A A 

mates of T and <p are necessary to achieve maximum output 
from the correlator/matched filter detector portion of the 
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Figure 11. Subprogram SGNLGN Flowchart 
receiver. That is, the maximum receiver sensitivity is 

A A 

obtained when ■ t = x and <p = <p as indicated in equations 

Pi 

(2.65) and (2.66). 

For the purpose of this study, the actual doppler 
shift is always set to zero in the transmit signal syn- 
thesizer, and the actual time delay due to the range or 
distance between array is set by the system geometry under 
consideration and' the reference speed of sound at the 
transmit array. Therefore, the estimate of doppler shift is 
set to zero in subprogram SGNLGN, and the estimate of time 
delay due to range is computed from the line-of-sight range 
between the center element of the transmit array and the 
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center element of the receive array, and the speed of sound 
c Q at the center element of the transmit array, or 



T 



V 




(3.10) 



where (x Q ,y o ,z o ) are the coordinates of the center element 
of the transmit array, and ( x r 'y r ' z r ) are the coordinates of 
the center element of the receive array. 



into the local processing waveform of equation (2.64) by 
a PPlying well-known properties of Fourier transforms to 
equation (2.17) to yield 



Note that the right hand side of equation (2.64) is just the 
time-sampled form of equation (3.11). 



(3.11) are identical to those used in equation (2.17) to 
generate the transmit signal in the ocean communication 
channel simulation computer program. Thus, the local 
processing waveform is identical to the transmit signal in 
functional form and total energy content, but is shifted in 
time and frequency by the estimates of range delay and 
doppler shift, respectively. 



The estimates of x and <j> are easily incorporated 



x (t-x) exp [ j 2 tt t ] = l c exp[j2TTgf ( t-x ) + j 2 tt 4> t ] (3.11) 

q=-K q ° 
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Subprogram SGNLGN also computes the local processing 
waveform signal energy E~ for later use in setting the 
decision threshold y. To compute the energy the magnitude- 
square of the complex Fourier coefficients are summed over 
all harmonics as indicated in equation (2.42) where c^ is 
the q-th harmonic of the local processing waveform. This 
sum is equivalent to the time-average power in the complex 
envelope of the local processing waveform. The energy can 
then be found by multiplying the time-average power by the 
fundamental period T q of the local processing waveform. 

3 . Subprogram AMPWGT 

The function of subprogram AMPWGT is to provide the 
real-valued, amplitude factors a m and b^ of the separable 
complex weights c m and d^ in equation (2.26). A flowchart 
of subprogram AMPWGT is shown in Figure 12. 

To generate the rectangular amplitude window, a m 
and b n are set equal to 1.0 for all elements (m,n) in the 
receive planar array. A separate subprogram to generate the 
amplitude weights facilitates generation of other forms of 
amplitude windows such as triangular, Hamming, Blackman, 
etc. However, only the rectangular window is used in this 
study . 

4 . Subprogram PHSWGT 

The function of subprogram PHSWGT is to generate the 
phase factors 0 m and cf> n of the separable complex weights 
c^ and d^ in equation (2.26) . A flowchart of subprogram 
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Figure 12. Subprogram AMPWGT Flowchart 

PHSWGT is shown in Figures 13a and 13b. PSHWGT can be 
programmed to compute phase factors that compensate, or 
remove the effects of, system geometry and deterministic 
medium wave propagation effects. PHSWGT can also introduce 
random noise in the phase factors for study purposes. 

The phase corrections for system geometry are com- 
puted using equations (2.27) and (2.28), where the direction 
cosines u and v are selected such that receive array beam 
pattern is aimed at the transmit array along the direct path 
from the receive array to the transmit array. That is, if 
n Q in equation (2.5) represents the direction of the direct 
path from the transmit array to the receive array, u and v 



71 




Figure 13a. Subprogram PHSWGT Flowchart 



72 





Figure 13b. Subprogram PHSWGT Flowchart 
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are set equal to -u and -v so that the receive array 

A 

beam pattern points in the -n direction. 

The deterministic medium phase correction factor 
is found by using equation (2.32) to compute the phase shift 
of the signal due to propagation through an inhomogeneous 
medium, and then negating the result. The random medium 
effect can also be computed from a closed form expression, 
but its use in generating phase correction factors is not 
the subject of this study, and will not be discussed. The 
total phase correction factor for the y-direction is 
obtained by adding the system geometry and deterministic 
correction factors. 

5 . Subprogram ARYPRO 

The array processor subprogram ARYPRO uses DFT and 
IDFT algorithms to: 

- generate the spectrum of the input electrical signal 
data at each element, 

- correct the phase of the spectral components to co- 
phase the signals at all array elements, and 

- inverse transform the signal spectrum at each array 
element to recover the cophased signal data. 

A block diagram of the subprogram ARYPRO is shown in Figures 

14a and 14b. 

The output of the array processor is a time-sampled, 
complex envelope signal representing the sum over all array 
elements of the signals at each element'. The defining 
equation for the DFT is used instead of a fast Fourier 
transform algorithm because the maximum number of time 
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Figure 14a. Subprogram ARYPRO Flowchart 
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Figure 14b. Subprogram ARYPRO Flowchart 
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samples permitted (33) is insufficient to achieve a measurable 
improvement in speed of execution [Ref. 6:pp. 151-152]. 
Additionally, use of the DFT equation relates the structure 
of the program directly to previous work by Ziomek [Ref. 2] 
and Vos [Ref. 3] . 

6 . Subprogram AWGN 

The function of subprogram AWGN is to generate one 
sample of an uncorrelated, Gaussian process with arbitrary 
mean and variance. A flowchart of AWGN is shown in Figure 15. 



Subroutine Library (IMSL) FORTRAN pseudorandom number 

generator routine GGNQF. GGNQF is a function subprogram 

that returns one zero mean, unit variance, Gaussian, or 

N(0,1) , pseudorandom number with each call to the function 

subprogram. The zero mean, unit variance pseudorandom 

number x is then scaled with the desired mean y and the 
n x 

desired standard deviation o x using the relation 



where the standard deviation is computed from the variance 
as 



AWGN is based on the International Mathematical 



X 




(3.12) 




(3.13) 



The desired mean and standard deviation are passed 



to AWGN as arguments from the main program. The mean is 




Figure 15. Subprogram AWGN Flowchart 



always set to zero for the purposes of this study and the 
standard deviation is found from equation (3.13) where the 
variance represents the desired power level in either the I 
or Q channel noise signal. 

The function subprogram GGNQF internally recomputes 
a new seed value for subsequent calls from a seed provided 
on the first call to the function. Since the first seed 
value is set as a parameter in the top level program, and 
is passed as an argument through AWGN to GGNQF, the pseudo- 
random sequence always follows the same pattern each time 



78 



the receiver simulalion is run. Usiny the same pseudo- 
random sequence for each simulation run is essential if 
comparisons between different runs are to be made from the 
receiver simulation output plots. 

7 . Subprogram INTGRT 

The function of subprogram INTGRT is to compute, 
using equation (2.68), the approximation of the correlation 
integral given in equation (2.67). A flowchart of subprogram 
INTGRT is shown in Figure 16. The complex product of the 




Figure 16. Subprogram INTGRT Flowchart 



total signal and noise, and the local processing waveform 
is computed in the top level program. The resulting 
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complex- valued sequence xs passed to INTGRT which computes 
the correlator output ~t by using a trapezoidal integration 
algorithm to approximate the integral of equation (2.67). 

A separate subprogram allows easy implementation of other 
numerical integration algorithms, if desired. 

8 . Subprogram WRITBL 

The function of subprogram WRITBL is to simply pro- 
vide tabular output of the simulation parameters and selected 
information generated during the execution of the program. 

A flowchart of subprogram WRITBL is shown in Figure 17. 




Figure 17. Subprogram WRITBL Flowchart 



Besides the data read into the COMMON storage blocks 
by subprogram READY, the information output to the data file 
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in table form includes the array element input SNR at which 
the Pd is estimated, the estimate of Pd, an estimate of Pfa 
obtained using equation (2.82), the SNR at the output of the 
array processor, the array gain and the SNR at the output of 
the magnitude-square operation. The tabular output is 
intended primarily to provide a convenient means of testing 
program modifications. It is not intended to be the primary 
simulation output, and will not be discussed further. 

9 . Subprogram PDPLOT 

The function of subprogram PDPLOT is to convert the 
numeric data generated by the receiver simulation into 
graphic form. A flowchart of the subprogram PDPLOT is shown 
in Figure 18. The plots generated from the numeric data 
are obtained through the use of standard DISSPLA graphics 
library subroutines. 

The primary output of the receiver simulation is a 
plot of the estimate of Pd versus the array element input 
SNR with a fixed value of Pfa as a parameter. PDPLOT also 
computes a theoretical value of Pd using the equation (2.82) 
where the SNR at the output of the magnitude-square operation 
is related to the array element input SNR through equation 
(2.91) . 

The curve representing the estimated Pd is plotted 
through the estimated Pd data using the DISSPLA least squares, 
cubic spline, curve fitting plot routine SMOOTH. The value 
of Pd obtained from equation (2.82) is plotted using the 
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Figure 18. Subprogram PDPLOT Flowchart 
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DISSPLA cubic spline, interpolating polynomial routine 
SPLINE, and appears as a dashed line on all plots generated 
by the receiver simulation. That is, the curve for the 
calculated value of Pd is an interpolating polynomial that 
passes through the calculated Pd data points, while the 
curve for the estimated value of Pd is allowed to pass 
within some arbitrary, small offset of the estimated Pd 
data point. For the purpose of this study, the maximum 
offset allowed was 10 percent of the data point value. 

This value of offset provided a reasonable compromise be- 
tween obtaining a relatively smooth fit of the plotted curve 
without diverging excessively from the data points. Both 
the estimated and calculated Pd data points are plotted at 
the same array element input SNR. 

B. MODEL VERIFICATION 

Verification of the receiver model involved two tasks, 
characterization of the pseudorandom number generator noise 
source, and a comparison of the data generated by the re- 
ceiver simulation program with the results predicted by 
the theoretical development discussed in Section II. Why 
simulate a test case for which a theoretical model exists? 
The main reason is to gain some confidence in the results 
generated by the simulation when the inputs are such that a 
theoretical model does not exist or is mathematically 
intractable . 
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1. Characterization of the Noise Source 



The capability to generate noise samples consistent 
with the noise description used in developing the receiver 
equations is central to verification and usefulness of the 
computer simulation. Thus, the first step is to ascertain 
the statistical properties of the noise source used in the 
computer program. Recall that the assumption required in 
developing the relations between Pd, Pfa, decision threshold 
and SNR of the magnitude-square of the correlator output 
involved the use of zero mean, Gaussian noise which is uncorre- 
lated in both spatial and temporal coordinates. The noise 
source was tested by taking a sequence of noise samples for 
both I and Q channels in manner identical to that used in 
the receiver simulation. That is, the subprogram AWGN 
was embedded in a test program to ascertain the statistical 
properties of the complex envelope of the noise signal. The 
following tests were applied to the I and Q channel sample 
sequences to verify agreement with the noise process 
assumptions : 

- sample mean 

- sample variance 

- histogram 

- autocovariance, and 

- estimated power spectral density. 

The test algorithms were obtained from standard IMSL 
procedures . 
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Figures 19 and 20 show the results of the sample 
mean, variance and histogram calculations for the noise 
processes in the I and Q channels, respectively. The sample 
mean and variance calculations approximate the zero mean, 
unit variance assumptions quite well. The histogram 
calculations indicate the distribution of the samples is 
roughly Gaussian. The smooth curve plotted through the 
histogram represents the exact Gaussian distribution for the 
size of the sample window used. The greatest deviation from 
Gaussian appears near the mean value which unfortunately is 
where most of the sample values lie. Thus, some difficulty 
in achieving a perfect correlation between simulated per- 
formance and theoretical results was anticipated, and in 
fact some deviation from theory at large values of Pd did 
occur . 

The space-time correlation properties of the sample 
sequence was measured by computing the autocovariance func- 
tion of 2000 samples over a total time (or space) displace- 
ment (lag) of 1000 samples. The autocovariance was computed 
using the IMSL routine FTAUTO . The correlation of adjacent 
samples, whether one or two or 1000 samples apart, was 
found to b.e remarkably small. The autocovariance function 
of the I and Q channel noise is shown in Figures 21 and 22, 
respectively. Because of the manner in which the samples 
are drawn in the simulation, the distinction of a sample 
being assigned to a particular channel at a particular 
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Figure 19 . Histogram of I-Channel Noise 
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Figure 21. Autocovariance Function of I-Channel Noise 
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Figure 22. Autocovariance Function of Q-Channel Noise 



element at a particular time is immaterial. The auto- 
covariance of the samples indicates that noise samples can 
be generated which are independent of sample time, array 
element location, and receiver channel. 

The estimated PSD was computed primarily to rein- 
force the results obtained from the autocovariance compu- 
tation. The IMSL routine, FTFPS, used to compute the 
PSD of the sample sequence implements an algorithm similar 
to that due to Welch [Ref. 7:pp. 553-554] in which the total 
sample record is partitioned into contiguous subrecords. 

Each subrecord is amplitude weighted with a triangule window 
function. Then, a periodigram of each amplitude weighted 
sample subrecord is computed. The resulting subrecord 
periodigrams are averaged over all subrecords in an effort 
to reduce the variance of the estimate of the PSD. The 
estimate of PSD for the I and Q channel noise is shown in 
Figures 23 and 24, respectively. The estimated PSD is 
approximately flat for both the I and Q channels over the 
sampled frequency range of 0 to n radians, or one sample 
frequency period. This tends to support the results obtained 
from the autocovariance computation that the noise samples 
are indeed uncorrelated. 

2 . Verification of the Output Data 

The primary output data from the receiver simulation 
is a plot of an estimate of Pd for various values of input 
SNR at a given Pfa. As shown in equation (2.82), a closed 
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ONE PERIOD OF SRMPLED SPECTRUM 
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form relation exists between the Pd, Pfa and correlator 
output SNR. Furthermore, the correlator output SNR can be 
related to the array input SNR when the signals at all array 
elements are precisely cophased by equation (2.91). 
Therefore, a plot of the theoretical performance computed 
using equations (2.82) and (2.91) can be superimposed on 
the plot of the receiver simulated performance. The curve 
representing theoretical performance provides verification 
of the simulated performance when input signal parameters 
and phase weighting allow precise signal cophasing to 
occur. The computed performance curve also provides a 
baseline to evaluate simulated performance when the input 
noise or receiver operating conditions differ from the 
assumptions used in formulating the receiver model. In 
all the plots, the theoretical performance is shown as a 
dashed curve. 
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IV. PRESENTATION OF RESULTS 



This section presents the results of three receiver 
simulation case studies and analyzes the simulation output 
in each case. The three cases considered are: 

- Case TEST, in which the output electrical signal data 
at each array element is produced by a test signal 
generator computer program, and is free of communi- 
cation channel effects such as attenuation, phase 
shifts due to wave front refraction, time delay due to 
range separation, and the transmit array beam pattern. 

- Case HMG1 , in which the output electrical signal data 
at each array element is produced by the ocean 
communication channel simulation computer program, 
and includes the effects of path attenuation, phase 
shifts due to system geometry, time delay due to range 
separation, and the transmit array beam pattern. 

- Case INHMG1 , in which the output electrical signal 
data at each array element is produced by the ocean 
communication channel simulation computer program, 
and includes not only the effects given in case HMG1 , 
but also contains phase shifts due to refraction of 
the wavefront along the propagation path. 

The results of these case studies will be analyzed by 
providing a brief description of the transmit signal 
waveforms used to generate the array element output elec- 
trical signal data, by listing the system parameters that 
distinguish the test cases, and finally, by interpreting 
the plots of Pd versus array element input SNR for each 
case study. 



A . TRANSMIT WAVEFORMS 

The analysis and theoretical development of the receiver 
simulation in Section II makes no assumption regarding the 
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functional form of the transmit signal. The receiver 
should perform equally well regardless of the type of trans- 
mit signal provided that the correlator/matched filter uses 
a replica of the transmit signal as the local processing 
waveform. To verify this hypothesis, and thereby test the 
integrity of the receiver simulation and the validity of the 
assumptions used in generating the model, more than one type 
of transmit signal was used to test receiver simulation 
performance. The use' of multiple transmit signal waveforms 
also provides a broader base from which to draw conclusions 
regarding the concept of model-based signal processing. For 
the purpose of this research, two types of transmit wave- 
forms were used in producing the Pd versus array element 
input SNR plots in each of the three case studies. 

A rectangular envelope CW pulse and a rectangular 
envelope LFM pulse were selected as transmit signals. These 
two particular forms were chosen for several reasons. First, 
the CW pulse and LFM pulse waveforms were considered to be 
typical transmit signals used in SONAR and acoustic signal- 
ing systems. Second, the CW pulse provides a simple case 
of amplitude modulation while the LFM pulse provides an 
example of an amplitude and angle modulated waveform with 
considerably different spectral characteristics. Finally, 
either signal can be readily synthesized from a finite, 
complex, Fourier series whose coefficients c^ may be calculated 
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from a closed form expression for the frequency spectrum 
of the time domain pulse characteristics. 

The rectangular envelope CW pulse is discussed first. 

The pulse repetition frequency is the same as the funda- 
mental frequency f of the finite frequency spectrum from 
which the pulse is synthesized. The pulse duty cycle was 
arbitrarily taken to be 0.5 yielding a pulse width of T q /2 
seconds where T q = l/f Q . The complex, Fourier series 
coefficients used to synthesize the complex envelope of the 
CW pulse are obtained from a closed form expression for the 
complex-valued continuous spectrum. To obtain the Fourier 
coefficients, the closed form expression for the continuous 
spectrum is evaluated at the discrete frequencies qf Q , and 
the resulting complex value is divided by the fundamental 
pulse period T . The index q takes on integer values 
between -K to K where K is the maximum number of harmonics 
in the finite Fourier series used to synthesize the CW 
pulse. The continuous spectrum of the CW pulse is the 
familiar sin(x)/x form obtained by taking the Fourier trans- 
form of the rectangular pulse shape where the pulse width 
is one-half the pulse period. The following specific 
transmit signal parameters were used in all CW pulse 
simulations : 

A = 40.0 

D = 0.5 

f 

o 



- Amplitude, 

- Duty Cycle, 

- Fundamental Frequency, 
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200 Hz 



- Number of harmonics, 

- Harmonic values. 



KMAX = 5 



c 

o 




= 20.00000 


exp [ j 0 ° ] 


C -1 


= C 1 


= 12.33240 


exp [ jO ° ] 


C -2 


= C 2 


= 0.000000 


exp [ j0°] 


C -3 


= C 3 


= 4.244132 


exp [ jl80°] 


C -4 


= C 4 


= 0.000000 


exp [ j0°] 


c -5 


= C 5 


= 2.546479 


exp [ j 0 ° ] 



The LFM pulse, or linear frequency modulated pulsed 
carrier waveform, is discussed next. The complex Fourier 
coefficients used to synthesize the LFM pulse are found 
using a procedure similar to that used for the CW pulse 
except the closed form expression for the complex-valued 
continuous spectrum of the LFM pulse was found by using the 
method of stationary phase. Officer [Ref. 8:pp. 67-68] 
describes the method of stationary phase as does Papoulis 
[Ref. 9:pp. 267-273] who also provides a complete descrip- 
tion of the LFM waveform. The following transmit signal 
parameters are used in all LFM pulse simulations: 



- Amplitude, 

- Duty Cycle, 

- Phase Deviation Constant 

- Fundamental Frequency, 

- Number of harmonics, 

- Harmonic values, 



A =40.0 

D =0.8 

B = 2356.2 radians/volt 

f = 10 Hz 

o 

KMAX = 3 

c = 14.60593 exp [ j 4 5 ° ] 

o 

c_ 1 = c 1 = 14.60593 expl j21°] 
c_ 2 = c 2 = 14.60593 exp[ j309° ] 
c_ 3 = c 3 = 14.60593 exp [ j 1 8 9 0 ] 
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It should be emphasized that these complex, Fourier 
series coefficients which are used in the ocean communica- 
tion channel simulation computer program to synthesize the 
complex envelope CW and LFM pulse transmit waveforms are 
also used in subprogram SGNLGN to produce the local processing 
waveform for the correlator/matched filter detector section 
of the receiver simulation computer program RCVR. 

B . CASE TEST 

The input signal data for case TEST is obtained from a 
separate computer program written to synthesize array 
element output electrical test signals from a finite Fourier 
series expansion of the test signal waveform. The technique 
is similar to that used in generating the receiver simula- 
tion local processing waveform in subprogram SGNLGN. How- 
ever, appropriate phase shifts are added to the Fourier 
coefficients to vary the direction of arrival of the plane 
wave incident on the receive array, and to incorporate the 
time delays due to spatial separation, d^ and d^, of the 
array elements. Furthermore, the time delay due to range 
and the doppler shift of the test signals are both set to 
zero to prevent signal loss in the correlator/matched filter 
due to a mismatch between the total array output signal and 
the local processing waveform. In this way, array element 
output electrical signals are obtained that exclude effects 
due to propagation of the sound wave through the ocean 
medium, or spatial modulation of the signals due to the 
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transmit array beam pattern. Since the medium effects are 
excluded from the case TEST array element signals, the 
carrier frequency f is set to zero, and the array element 
spacing is increased to keep the interelement spacing d^ 
and dy equal to one-half the shortest wavelength in the 
transmit signal. Note that this produces different inter- 
element spacings for the CW and LFM pulse situations. 

Signals generated in this fashion allow the receiver 
output to be analyzed separately from the communication 
channel. Essentially, it is equivalent to having a signal 
generator that provides a test signal to measure receiver 
performance. The results obtained from case TEST, when 
combined with the plot of theoretically predicted performance, 
provides a baseline with which to analyze the receiver 
output when the ocean communication channel simulation data 
is processed. 

The system parameters defining case TEST are as follows: 

- Array Parameters 

Number of array elements, M = N =5 

Array Element Spacing, CW d^ = d^ = 0.7375 meters 

LFM d = d =24.58 meters 
x y 



c q = 1475 meters/second 
= 0.0 seconds 



- Medium Parameters 
Speed of Sound, 

Actual time delay 
Actual doppler shift ^A = Hertz 

Four plots representing receiver performance were 
generated for case TEST. Phase weighting is done to 
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compensate for the direction of arrival of an acoustic plane 
wave incident upon the array elements. The direction of 
arrival can be set to any arbitrary direction by the test 
signal generator computer program. However, the direction 
was chosen to be the same as that due to the system geometry 
used in the case studies HMG1 and INHMG1 . Phase weighting 
for geometry is indicated by the state of the logical 
variable STEER. The estimates of time delay due to range 
and doppler shift were both set to zero in subprogram SGNLGN 
to maximize the output of the correlator/matched filter 
detector. Two of the plots. Figures 25 and 26, show receiver 
performance for a rectangular envelope C W pulse. Figure 25 
was produced using a Pfa of 0.1, and Figure 26 shows the 
effect of a Pfa of 0.01. The remaining two plots generated 
for case TEST, Figures 27 and 28, show receiver performance 
for a rectangular envelope LFM pulse at a Pfa of 0.1 and 
0.01, respectively. Since the input electrical signal data 
does not incorporate channel effects, phase weighting for 
deterministic and random medium effects was not done for 
case TEST. 

Several observations regarding model performance can 
be made by analyzing the plots. Generally for all plots, the 
measured performance of the simulation (solid curve) follows 
the predicted theoretical performance (dashed curve) , but 
the agreement is not exact. That is, the estimate of Pd 
follows the same trend as the predicted Pd, but the 
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SIMULATED RECEIVER PERFORMANCE 
AS A FUNCTION OF 
SNR AT A SINGLE ARRAY ELEMENT 
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transition from a low Pd to a Pd of 1 is more sensitive to 
changes in the array element input SNR. A comparison 
between the plots shows the Pd to be a function of the Pfa 
as predicted by equation (2.82), and for small values of 
SNR, the estimated Pd approaches the Pfa parameter used by 
the simulation in generating the plots. This effect is 
also predicted by equation (2.82) since at SNR values much 
less than one, the input signal is essentially off, and only 
the noise is producing signal detections. A comparison 
between the plots for the CW pulse and LFM pulse does show 
some minor difference in receiver performance, but the same 
general trends seem to hold regardless of the functional 
form of the transmit signal. 

The differences in the measured performance between the 
CW pulse and the LFM pulse simulations can be attributed in 
part to the difference in the maximum number of harmonics 
used to synthesize the transmit signals. As shown by 
equation (2.22), the number of time samples L is related to 
the maximum number of harmonics K used to synthesize the 
transmit signal. The difference in the number of time 
samples causes different subsequences of noise samples to be 
drawn from the noise generator subprogram AWGN. The 
different noise sample subsequences will cause small varia- 
tions in the relative frequency estimate of Pd. Thus, an 
exact comparison of the noise performance cannot be made 
across different types of transmit signals when the number 
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of harmonics used to create the transmit signal are not the 
same. However, the real significance of this effect is that 
even with different noise sample subsequences in the I and 
Q receiver channels at the array elements, the measured 
performance conforms closely to the performance predicted by 
theory under the very ideal assumptions of precise signal 
cophasing, and zero-mean, white, Gaussian input noise. 

For these reasons, it can be inferred that the basic 
receiver model is valid, but some difficulty may exist in 
generating precisely white, zero-mean, Gaussian noise within 
the simulation program. Indeed as indicated in Section 
III.B.l, the noise source used in the receiver simulation 
does not generate precisely zero-mean, white, Gaussian 
noise. Furthermore, the measured Pd is a relative frequency 
estimate of the actual receiver performance, and may also 
contribute some error in the simulation results. 

Considering these approximations to the assumptions used 
in the receiver model development of Section II, exact 
agreement should not be expected even in the case where the 
communication channel effects are excluded. However if 
only relative changes or effects are to be observed, moderate 
disagreement between the theory and simulated receiver 
performance can be negated by comparing only differences in 
the estimate of Pd for each of the case studies. A com- 
parison of the differences in the measured receiver perfor- 
mance will indicate if, and to what degree, the effects due 
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to the physics of acoustic wave propagation can be compen- 
sated for by model-based signal processing. The results of 
case TEST then can be viewed as a validation of the receiver 
model and the applied programming used to implement the 
model . 

C. CASE HMG1 

The input signal data for case HMGl and case INHMG1 is 
produced by the ocean communication channel simulation 
program due to Vos and Ziomek [Ref. 3]. For case HMGl, the 
ocean is considered to be a homogeneous medium with respect 
to the speed of sound propagation through the water. That 
is, the speed of sound is identical at the transmit and 
receive arrays, and at all points in between. In this model 
of the ocean medium, refraction or ray bending of the sound 
wave does not occur. Thus, only phase weighting for geometry 
is required, and only those results will be presented. The 
homogeneous case provides a needed baseline with which to 
judge the effects of model-based signal processing when 
inhomogeneous case data is studied. 

Since case HMGl includes the effects of the ocean com- 
munication channel, a carrier frequency is used to convey 
the modulation through the channel. In case HMGl and case 
INHMG1 the carrier frequency was set to 5.0 kHz. The array 
interelement spacing was adjusted accordingly to maintain 
the one-half wavelength spacing needed to prevent grating 
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lobes in the array beam patterns. Note that the interelement 
spacing was set using the highest frequency expected for 
either the CW pulse or the LFM pulse waveforms, and was 
kept constant for both the HMG1 and INHMG1 case studies. 

The highest frequency occurred for the CW pulse case, and was 
equal to 6000 Hz. The estimates of the time delay and 
Doppler shift in subprogram SGNLGN were set equal to the 
actual values to maximize the correlator output. 



The system parameters defining 
- Array Parameters 


case HMG1 are as follows: 


Number of array elements, 


M = N = 5 


Array Element Spacing, 
- Medium Parameters 


d = ^ 0.1229 meters 

x y 


Speed of Sound, 


c q = 1475 meters/second 


Actual line of sight 
time delay 


t a = 2.033898 seconds 


Actual Doppler shift 
- System Geometry (See Figure 1) 


<j> A = 0.0 Hertz 


Depth of Transmit Array 


y = 1000.0 meters 


Depth of Receive Array 


y = 2500.0 meters 
1 r 


Cross Range 


x -x = 500.0 meters 
r o 


Line of sight range 


|r-r 1 = 3000.0 meters 
1 o 1 


Eight plots were generated for 


case HMG1. The first four 


plots, Figures 29 through 32, were 


produced using the 



rectangular envelope CW pulse transmit signal. The remain- 
ing four plots. Figures 33 through 36 were produced by 
the rectangular envelope LFM pulse transmit signal. Figures 
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-40.0 -35.0 -30.0 -25.0 -20.0 -15.0 -10.0 -5.0 

INPUT SIGNAL-TO-NOISE RATIO (dB) 

Figure 29. Receiver Performance for case HMGl, CW Pulse, 

Pfa = 0.01, No Phase Weighting 
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29, 31, 33 and 35 show receiver performance when no phase 
weighting is done. These plots indicate that performance 
will be very poor when the output signals from the array 
elements are not cophased. Again, it is significant that 
the Pd in this case is approximately equal to the value of 
the Pfa parameter used in obtaining the plot. That is, 
the total array output signal must be very small because 
of destructive interference of the array element signals, 
and the only detections appear to be due to the noise. 

Figures 30, 32, 34 and 36 show performance when phase 
weighting is done to compensate for system geometry. These 
plots indicate receiver performance can be made to approxi- 
mate the theoretical predicted performance if the individual 
array element signals can be cophased to provide constructive 
interference at the output of the array processor. 

Figures 29, 30, 33 and 34 show receiver performance at 
a Pfa of 0.1. Figures 31, 32, 35 and 36 show the effect of 
decreasing the Pfa to 0.01. The effect of decreasing the 
Pfa on the simulated performance follows the theoretically 
predicted performance and provides further confirmation that 
the computer implementation of the model is accurate. 

The important conclusion from case HMG1 is dramatic 
improvement in receiver performance is possible provided 
that enough information about system geometry exists to 
generate phase weights that will properly cophase the array 
element output signals. 
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D. CASE INHMG1 



The input signal data for case INHMG1 is produced in 
the same manner as case HMG1. The transmit signal waveforms, 
system geometry, and transmit and receive arrays are identi- 
cal. The only difference between the two cases is the 
model of the ocean medium. For case INHMG1, the ocean is 
considered to be an inhomogeneous medium with respect to the 
speed of sound propagation through the water. The speed of 
sound is assumed to vary linearly with depth (or y-coordinate) 
only. That is, the speed of sound is not identical at the 
transmit and receive arrays. In this space-variant model 
of the ocean medium, refraction or ray bending of the sound 
wave does occur, and the model of the ocean medium due to 
Ziomek [Ref. 2] predicts what the phase shifts in array ele- 
ment output electrical signals will be. Using this knowledge 
of the physics of acoustic wave propagation, phase weights 
are generated that attempt to compensate for the refraction 
of the acoustic wave along the ray path, and the system 
geometry. To determine the effect of this model-based signal 
processing algorithm, the inhomogeneous data was also 
processed using phase weights that compensate just for sys- 
tem geometry. The difference in simulated receiver perfor- 
mance between these two situations will graphically show the 
effectiveness of the model-based signal processing approach. 

The system parameters defining case INHMGl are as 
follows : 
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- Array Parameters 

Number of array elements, 
Array Element Spacing, 

- Medium Parameters 

Speed of Sound, 

Gradient 

Actual line of sight 
time delay 



M = N = 5 



d = d = 0.1229 meters 
x y 



c q = 1475 meters/second 
q = .017 



x A = 2.033898 seconds 



<|> A = 0.0 Hertz 



y = 1000.0 meters 

J o 

y = 2500.0 meters 
1 r 

x -x = 500.0 meters 
r o 



= 3000.0 meters 



Actual Doppler shift 
- System Geometry (See Figure 1) 

Depth of Transmit Array 
Depth of Receiver Array 
Cross Range 

Line of sight range |r -r 

Twelve plots were generated for case INHMG1. The first 
six plots. Figures 37 through 42 show receiver performance 
for a rectangular envelope CW pulse. The remaining six 
plots. Figures 43 through 48, are the result of a rectangu- 
lar envelope LFM pulse transmit signal. Differences in 
receiver performance due to the form of the transmit signal 
are measurable, but not significant. 

Figures 37, 40, 43 and 46 were generated when no phase 
weighting is applied to the array output electrical signals. 
Again, these figures indicate that the array element output 
electrical signals must be cophased to provide any useful 
detection capability in the receiver. 
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SIMULATED RECEIVER PERFORMANCE 
AS A FUNCTION OF 
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SIMULATED RECEIVER PERFORMANCE 
AS A FUNCTION OF 
SNR AT A SINGLE ARRAY ELEMENT 
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Figures 38, 41, 44 and 47 show the receiver performance 
when only the system geometry is compensated for by the 
phase weights used in the array processor. These figures 
indicate that in the presence of the inhomogeneous ocean 
medium, phase weighting to compensate for system geometry 
will not achieve the receiver performance predicted by 
theory. That is, traditional beam steering is not suffi- 
cient for maximum receiver performance, and a large margin 
for improved performance exists when the effects of the 
ocean medium on the received signal can be predicted. Note 
also that the degradation in performance is consistent for 
the different transmit waveforms and values of Pfa. 

Figures 39, 42, 45 and 48 show the effect on receiver 
performance when phase weights are computed based on both 
the system geometry and the refraction of the acoustic wave 
in the ocean medium. The state of the logical variable 
DMEDIA indicates if the phase weights contain the correction 
for deterministic, inhomogeneous medium wave propagation 
effects. As expected, the simulated receiver performance 
conforms to the performance predicted by theory, and is 
nearly identical to that of case TEST when all medium 
effects have been corrected for. 

However, the really significant result is seen by com- 
paring Figures 38, 41, 44 and 47 with Figures 39, 42, 45 and 
48. The only difference between these plots is in use of 
phase weights that compensate for the phase shifts due to 
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propagation through an inhomogeneous ocean medium, and the 
model of the inhomogeneous ocean medium is based on the 
physics of acoustic wave propagation. If one assumes a 
marginal detection probability of 50% (Pd = 0.5), the 
improvement in receiver performance due to the application 
of a model-based signal processing approach for this particu- 
lar simulation geometry and ocean medium characteristic 
varies from 14 dB to 18 dB depending on the type of 
transmit waveform or value of Pfa. The conclusion is that 
the physics of wave propagation can be used to develop 
signal processing algorithms having significant impact on 
the performance of a receiver designed to process complex 
envelope, electrical signals from an array of point source 
elements . 
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V. CONCLUSIONS AND RECOMMENDATIONS 



The original objectives of this research have been accom- 
plished. The computer simulation of a correlator/matched 
filter receiver to process data from an array of electro- 
acoustic transducers has been developed and validated. The 
concept of model-based signal processing was applied to the 
development of an array signal processing algorithm, and was 
shown to have a marked impact on the capability of the 
receiver to detect the presence of a signal in zero-mean, 
additive, white Gaussian noise. 

The computer simulation of the receiver was validated by 
a direct comparison between the measured performance of the 
simulation and the performance predicted by theory when all 
array element output electrical signals are precisely 
cophased. The agreement between predicted and measured 
performance was close but not exact since the noise assump- 
tions made in developing the theoretical performance could 
not be precisely duplicated by the computer simulation. 
However, the close agreement with theoretical performance 
was considered to be a validation of the receiver simulation 
computer program. 

Test cases were designed and run that permitted relative 
changes in receiver performance to be measured. In this 
way, the effectiveness of the model-based signal processing 
approach could be quantified. For the simulation parameters 



134 



and transmit signals used in this study, compensating for 
the refraction of the acoustic field in the ocean medium 
was found to improve receiver performance over traditional 
beam steering by at least 14 dB when measured at the point 
of 50% probability of detection. 

Suggestions for further study of the receiver model 
include : 

- Obtain, or develop, a pseudorandom noise generator 
that more closely approximates the noise model 
assumptions implicit in the derivation of the 
theoretical receiver performance, and revalidate 
the computer simulation. 

- Simulate more realistic noise processes, such as 
colored noise, and determine the impact on receiver 
performance . 

- Try decision rules other than Neyman-Pearson , such 

as the minimum average probability of error criterion 
to measure probability of error performance for an 
underwater, acoustic data communication system. 

- Experiment with other functional forms for the 
transmit signal, and measure the effect on perfor- 
mance when errors in the estimate of time delay and 
Doppler shift exist. 
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